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O . Abstract 

The effect of molecule size (excluded volume) and the range of interaction on the surface 
^ tension, phase diagram and nucleation properties of a model globular protein is investigated us- 

ly^ ■ ing a combinations of Monte Cai^lo simulations and finite temperature classical Density Func 

2 ■ tional Theory calculations. We use a parametrized potential that can vary smoothly from the 

o : 

standard Lennard- Jones interaction characteristic of simple fluids, to the ten Wolde-Frenkel 
model for the effective interaction of globular proteins in solution. We find that the large ex- 

X 

' eluded volume characteristic of large macromolecules such as proteins is the dominant effect 

in determining the liquid-vapor surface tension and nucleation properties. The variation of the 
range of the potential only appears important in the case of small excluded volumes such as 
for simple fluids. The DFT calculations are then used to study homogeneous nucleation of the 
high-density phase from the low-density phase including the nucleation barriers, nucleation 
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pathways and the rate. It is found that the nucleation barriers are typically only a few kgT and 
that the nucleation rates substantially higher than would be predicted by Classical Nucleation 
Theory. 

March 31, 2010 

Introduction 

One of the most important problems in biophysics is the characterization of the structure of pro- 
teins. It is well-known that the main impediment to the determination of protein structure is the 
difficulty with which good quality protein crystals can be produced. This has led to a large body of 
work focused on the understanding the details of protein nucleation. In recent years, it has become 
apparent as a result of simulation, theoretical and experimental studies that nucleation in general, 
and protein nucleation in particular, is strongly affected by the presence of intermediate, metastable 
states.-"- This raises the possibility that by understanding the mechanism by which intermediate 
states affect nucleation, the quality of the final result can be better controlled. 

The practical importance of investigating the role of intermediate metastable states lies in the 
fact that the effective interactions of proteins in solution depend on their environment. Protein 
molecules interact via Coulombic forces mediated by the ionic solution in which they are dissolved. 
The effective interaction between two protein molecules therefore depends on the properties of the 
solution, particularly the salt used to create the solution and its pH. This is the reason for the well- 
known fact that some salts are more effective than others in precipitating protein crystallization 
(the Hofmeister effect).^ Detailed confirmation of the connection between the properties of the 
solution and the effective intra-protein interactions has come from both computer simulation-^ 
and from theoretical studies of the phase diagram of proteins in solution.- One goal of the present 
work is to investigate what aspect of the effective interactions is most relevant in controlling the 
nucleation of the metastable phase. 

Assuming that the effects of the solvent can, at first approximation, be entirely accounted for by 
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an effective interaction potential between protein molecules, the problem of protein nucleation can 
be viewed as analogous to the nucleation of a solid from a dilute gas. In this way, Wilson observed 
that favorable conditions for protein nucleation are correlated with the behavior of the osmotic 
virial coefficient.— Rosenbaum, et al. showed that the phase diagram of a large class of globular 
proteins can be mapped onto that of simple liquids with an interaction potential that depends on 
the ionic strength of the solvent.-^^^ One particular characteristic of proteins is that these models 
involve a hard-core repulsion and an attractive tail with the range of the attraction being quite 
small compared to the size of the hard core. It is known that as the range of the attraction becomes 
smaller, the critical point of the liquid-vapor transition is suppressed relative to the triple point 
until for sufficiently short-ranged potentials, the liquid- vapor transition becomes metastable with 
respect to the vapor-solid transition.- It is in this circumstance that the meta-stable liquid phase is 
thought to play an important role in protein nucleation. 

Despite the abundance of evidence for the role of the intermediate state from both simulation 
and experiment, there is still no really convincing theoretical description based on first-principles. 
Lutsko and Nicolis showed that transitions from vapor to metastable liquid to solid appeared ad- 
vantageous relative to the direct vapor-solid transition based on the bulk free energy surface,— 
but that work neglected the effect of interfaces and, in particular, surface tension. The goal of 
the present work is to take a step towards filling in this gap by characterizing the liquid- vapor in- 
terfacial surface tension and the liquid-vapor transition in a model globular protein consisting of 
molecules interacting with the ten Wolde-Frenkel potential.- However, because of the metastable 
nature of the transition, such a characterization based on simulation is difficult. We found it impos- 
sible to stabilize the liquid- vapor interface for the model protein due to the strong tendency towards 
crystallization. We have therefore had to use an indirect approach consisting of a combination of 
theory and simulation. First, the ten Wolde-Frenkel model potential for globular proteins is gener- 
alized so that it depends on three independent parameters that allow it to be deformed continuously 
from a Lennard- Jones potential, i.e. a simple fluid, to the hard-core-i-tail ten Wolde-Frenkel (tWF) 
potential. We have performed simulations covering part of the range from simple fluid to protein 
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and compare this to Density Functional Theory calculations to show that the DFT remains quan- 
titatively accurate as the hard-core radius is increased from zero, in the LJ potential, to a typical 
value in the tWF potential and as the range of attraction decreases. Whereas simulation becomes 
infeasible when the range becomes too short, the calculations are easily performed and the preced- 
ing agreement gives some confidence in the result. The DFT is then used to calculate the nucleation 
barrier using recently developed energy surface techniques. We find, somewhat surprisingly, that 
the increase in the hard-core radius is much more important than the decrease in the range of the 
potential and is mostly responsible for a dramatic drop in surface tension. This observation is 
particularly relevant in the case of proteins since the properties of the effective interaction can be 
varied by changing controllable parameters, such as the PH of the solution, so that a primary goal 
is to determine those conditions most favorable to homogeneous nucleation. Our results indicate 
that, so far as the low-density /high-density part of the transition is concerned, varying the range of 
the potential has little effect. 

In the next Section, our model potential is defined and the Monte Carlo simulations described 
and the DFT calculations are also sketched. In Section III, we describe the comparison between 
theory and simulation as the potential is varied from Lennard- Jones towards the model protein 
interaction and show quantitative agreement between theory and simulation. We then present the 
DFT results for the surface tension and nucleation barrier and nucleation rate as calculated from 
DFT and compare to classical nucleation theory. The paper concludes with a discussion of our 
results. 

Theoretical Methods 

Simulation 

Because of the complexity of the constituent particles and of the solvent-induced interaction glob- 
ular proteins are often modeled, in a first approximation, by an effective interaction potential. In 
their study of the phase behavior of globular proteins in solution ten Wolde and Frenkel proposed 
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the following effective potential- 

This is a Lennard- Jones potential modified to capture the essential features that characterize the 
interaction of proteins in solution. The potential includes a hard-sphere repulsion at r = a that 
accounts for the excluded volume, i.e. the size, of the protein molecule. The parameter a controls 
the range of solvent induced effective attraction in that the minimum of the potential well occurs 



/ 1/3 

at r = ay 1 + (^) while the minimum of the potential is — £, independent of the value of 
the range. The range of the potential decreases when a increases and for a = 50 this model 
potential reproduces the phase behavior of globular protein solutions and with in particular the well 
known fact that for sufficiently short range attractions the liquid phase becomes thermodynamically 
unstable,2j5. as illustrated in [[figure] [1][]1[ For the potential (??) this behavior is expected to be 
found when a ~ 10 in accordance with the generally accepted criteria according to which the range 
of the attraction should be ~ 25% of the range of the repulsive part to have a metastable liquid-gas 
coexistence.- 

Metastable systems are notoriously difficult to study numerically because the liquid state spon- 
taneously decays to the more stable crystal phase. In the present situation, the value a = 50 puts the 
liquid-gas coexistence curve deep inside the unstable region which makes the liquid phase difficult 
to maintain. A possible solution to this problem is to use constraints that restrict the configurational 
space of the system thereby preventing the transition to the solid phase. There are different ways 
to constrain the system.— In the restricted Monte Carlo method— the density is constrained to be 
below a limited value which thus suppresses the dense phase. While this approach is well suited to 
the study of the supersaturated vapor phase it is not appropriate for studying supercooled liquids 
because the densities of the liquid and solid phases are too close. Another possibility would be 
to limit the number of neighbors per particle. This approach has been applied with success to a 
Lennard- Jones potential. However there is a body of evidence that besides its effect on the stability 
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Figure 1: Calculated phase diagrams as a function of a for 5 = 1 showing that the liquid is 
metastable a > 10. From Lutsko and Nicolis.- 
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of liquids, the range of the attractive potential impacts also on the structure of the liquid. This fact 
could be of special importance in the present case because the structure of the liquid is atypical due 
to the large excluded volume which limits the number of interactions per particle. More generally 
the shortcomings of restricted ensemble methods is the elimination of many configurations. 

Li this work no constraints were imposed on the system. Stable liquid-gas equilibrium were 
obtained starting from a stable two-phase state and modifying step by step the temperature and the 
parameter a until the desired conditions are met. Unfortunately for the desired case of a = 50 this 
was not possible because there is no stable configuration to start from. We therefore generalized 
the potential as 



with a hard core at r < 5o, so that, by varying the dimensionless parameters a and 5, it is possible 
to smoothly go from a simple fluid, described by the Lennard- Jones interaction with a = 1 and 5 = 
to the model protein interaction with 6 = 1 and a = 50. This allows to reach a metastable liquid- 
gas coexistence from a stable Lennard- Jones system by varying in small steps the temperature and 
the parameters oc and 5. 

We simulated an ensemble ofN= 1885 particles interacting via the potential (??) with a stan- 
dard Metropolis Monte-Carlo algorithm (MC-NVT).— The potential is truncated at rc = 2.8 but 
not shifted. The particles are contained in a volume V with dimensions Lj^ = Ly = 9a, and 

= 108a and periodic boundary conditions are imposed in all directions. To avoid the prob- 
lem of stability of the liquid phase, the system is equilibrated during 5x10^ Monte-Carlo cycles 
(one cycle = A'^ updates) at given (T, a, 5) starting from a stable initial configuration at (T + AT, 
a + Aa, 5 +A5). The configuration of the system consists in a liquid slab of thickness Az — 21 o 
surrounded along the z-direction by two gas slabs. In this way the stability of the liquid phase was 
maintained during the 10^ cycles used to measure the density profile and the surface tension. 

Surface tension can be measured by different methods. In a recent paper we adopted the Ben- 
nett method as this method appears to be precise.— However we found it difficult to implement 
here for the following reason. In the Bennett's method the calculation of the surface tension follows 




(2) 



7 



from the definition 

where F is the free energy and A = x is the area of each liquid-vapor interface. In its imple- 
mentation the method requires that one compares the energies (/ = 0, 1) of two configurations 
with different liquid- vapor interface areas Aq and A\. The configuration with interface area Ai is 
obtained from the previous one by rescaling the positions of the particles—:— 

x' = x(Ai/Ao)5, 
y = 3;(Ai/Ao)i 

z' - z(Ao/Ai). (4) 

Because this perturbation is done at constant volume, the system is expanded along one direction 
and compressed along the transverse direction. Due to the isotropy of the fluid phase this trans- 
formation has only a negligible effect on the energy of homogeneous systems because the energy 
change created by the displacement along the compression direction is compensated by an energy 
change along the expansion direction. In the presence of an interface the symmetry is lost and there 
is a net energy difference between the two configurations which is localized at the liquid-vapor in- 
terface. While the perturbation ^ is usually an efficient way to probe the interfacial free-energy 
we found it to be unsuitable here as the surface tension of the protein model is particularly small so 
that large perturbations are necessary to have a measurable free energy difference. Unfortunately 
large perturbations are inefficient when applied to systems whose potential contains an excluded 
volume because the compression step creates configurations with overlaps between particles. This 
not only gives infinite energy variation but these variations are often located inside the liquid phase 
and not at the interface. To solve this problem we implemented the interface wandering approach 
which allows a precise evaluation free-energy difference with small area perturbations.— In this 
method the system is perturbed at each Monte Carlo cycle according to (01). Contrary to the Ben- 
nett's method where the new configuration is only tested, in the interface wandering method the 
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new area, chosen at random in an interval [Amm j^max], is effectively accepted with probability 



^accepted = min ( 1 , e 



) 



(5) 



As a result of this acceptance ratio, the area of the interface evolves and is distributed between A, 
and Amax according to a distribution /, which is related to the free energy 



Because in this method many areas are sampled, instead of one in the Bennett's method, a curve 
fitting can be performed on the interface area distribution function which allows a precise determi- 
nation of the surface tension. 

DFT calculations 

The properties of inhomogeneous systems were calculated using Density Functional Theory. Ac- 
cording to DFT, the free energy in the grand canonical ensemble, the grand potential Q., is a func- 
tional of the local density, p (r) , and the applied external field, (j) (r) , of the form 



where /i is the chemical potential, F and Q. both depend on temperature and the functional F [p] 
does not depend on the field.— >^ The equilibrium density distribution is determined by minimizing 
Q. [p]. In a uniform system, the local density is a constant, p (r) = p and F[p] F (p) is the 
Helmholtz free energy. 

Li our calculations, F [p] is approximated using the Modified-Core van der Waals DFT model.— 
This is a generalization of the simplest hard-sphere plus mean-field tail model which gives quanti- 
tatively accurate descriptions of fluid structure, surface tension— and nucleation properties.— 
Here, we require quantitative accuracy in the DFT calculations as we wish to make direct compar- 



F = -ksTlnf 



(6) 




(7) 
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ison to simulation. In this model, the free energy is expressed as a sum of three contributions: 



F [P] = Fhs ( [p] ; dHs) + Fcore ( [p] ; dns) + Ftari ( [p] ; dHs) ■ (8) 

The first term on the right is the hard-sphere contribution which is treated using the Fundamental 
Measure Theory functional. The second term is the "core correction" which is similar in form to 
the hard-sphere term, but which modifies the hard-sphere contribution so that the free energy in 
the bulk phase reproduces a given equation of state.— The last term is the tail contribution and has 
the simple mean field form 



Ftaii{[p];dHs) = ^J ®{ri2-dHs)p{ri)p{r2)v{ri2)dridr2. 



(9) 



The reference hard-sphere diameter, dns, is calculated using the Barker-Henderson— expres- 
sion, 

dHS = £ {l-e-^''^'-^)dr (10) 

where ro is the distance at which the potential vanishes, V (ro) = 0. 

The model requires as input the bulk equation of state. For this, we use the Barker-Henderson 
first-order perturbation theory,— 

F{p)=FHs{p\dHs) + ^P^ j &{ri2-ro)v{ri2)gHs{n2;dHs)dridr2, (11) 
where gns {ri2',dHs) is the hard-sphere pair distribution function in the fluid phase. 

Results 

Surface Tension 



[figure] [2] [] 2 shows the phase diagrams obtained from the simulations as described above and as 
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Figure 2: (Color online) The phase diagrams for a = 1, panel (a), and for a = 5, panel (b), as 
determined by simulation and theory for different values of the hard-core radius, R = 6o. The 
curves, from top to bottom, are for 5 = 0,0.2,0.4,0.6,0.8 and 1.0, respectively. 

predicted using the thermodynamic perturbation theory for the cases of a = 1 and a = 5. Away 
from the critical point, the agreement between theory and simulation is satisfactory over a wide 
range of hard-core radii. As is usual for a mean-field theory, the agreement near the critical point 
is not expected to be very good and this region has therefore not been studied. The critical density 
and temperature for each system was estimated following the procedure in Ref. and are given in 



[table] [ 1 ] [] 1 1 and [[table] [2] []2| The coexistence curves are shown in [figure] [3] [] 3 with the density 
and temperature scaled to the critical density and critical temperature respectively. Although the 
generalized potential involves two length scales, o controlling the position of the minimum and 
the hard-sphere radius 5o, the coexistence curves show typical corresponding- states as has been 
seen for other multi-length scale potentials.— 

A comparison between the surface tension as determined by simulation and from DFT calcu- 



lations is shown in [figure] [4] [] 4 Increasing the size of the molecule, i.e. the excluded volume 



parameter 5, at fixed a leads to lower surface tension: since the surface tension scales, roughly, 

11 



with the critical temperature, see|[table][l][]l[ this trend is attributable to the decrease of the crit- 



ical temperature with increasing d which, in turn, is due to the increase with 5 of the range of 
repulsion (which extends from r = to r = do + \J 8^ + (^) ^^^) compared to the range of the 
attractive part of the potential (and taking into account that the depth of the potential is fixed). On 
the other hand, for fixed 5, increasing a has little effect for 5 = 1 but leads to a significant increase 
in the surface tension for 5 = 0. In the latter case, the potential can be written as a function of the 
single parameter <5 j a^l^ so that one expects the surface tension in this case to vary as 7 ~ a}l^ . 
(Note that in fact this scaling is somewhat spoiled by the fact that the cutoff used in the simulations 
was not scaled in the same way, so this should only be taken as an explanation for the general 
trends.) In the former case of 5 = 1, the effect of changing a is mitigated by the relatively large 
(and fixed) effect of the excluded volume and, of course, the fixed depth of the potential minimum. 

Since the surface tension goes to zero at the critical point and, as discussed above, the pertur- 
bative equation of state is least accurate near the critical point, it is not surprising that the relative 
accuracy of the theoretical calculations decreases as one approaches the critical point. Away from 
this region, however, the theoretical calculations are in good, nearly quantitative agreement with 
the simulations. Based on this comparison, the extension of the calculations to higher values of a 
seems justified. 



[table][l][]l| and |[table][2][]2| also show the result of a fit to the expected form 7 = 70 (1 — 
T /TcY'^^. Fitting both 70 and Tc and comparing to the critical temperature extracted from the 
coexistence data gives a consistency check on the analysis. In general, the temperatures derived 
by both methods are in reasonable agreement. The errors in the critical density are too large to 
permit us to make a definitive statement concerning whether or not the surface tensions obey a law 
of corresponding states. 

Nucleation of globular proteins 

Globular proteins are modeled using a = 50 and 5 = 1 so that the liquid is metastable as shown 



in [figure][l][]l As discussed above, this makes simulation extremely difficult so that we only 
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Figure 3: (Color online) The same as [figure] [2] [] 2 with the density and temperature scaled to the 



critical density and temperature for each potential. 

Table 1: Fits to surface tension: y{T) = 7o71(l — T /Tc)^-'^^ for a = l.The last two columns give 
the critical properties estimated from the coexistence data. 
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for a = 5. 



5. The last two columns give 
the critical properties estimated from the coexistence data. 
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Figure 4: (Color online) The surface tension as a function of temperature for a = I, panel (a), 
and for a = 5, panel (b), as determined by simulation and theory for different values of the hard- 
core radius, R = do. The curves, from top to bottom, are for d = 0,0.2,0.4,0.6,0.8 and 1.0, 
respectively. 
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Figure 5: (Color online) The surface tension as a function of temperature for a = 50 as calculated 
using DFT. The curves, from top to bottom, are for 6 = 0,0.2,0.4,0.6,0.8 and 1.0, respectively. 

present the results of calculations using DFT. Indeed, the fact that the metastable phase immedi- 
ately tends to nucleate the solid phase in simulation is one of the main reasons for focusing on 
the nucleation of the metastable phase as this appears to be the rate-limiting step. The agreement 
found above between DFT and simulation for lower values of a suggests that these results should 



be reasonable quantitative estimates, [figure] [5] [] 5 shows the surface tension as a function of tem 



perature for different values of the hard core radius. The trends noted earlier are repeated with the 
surface tension for 5 = being approximately 50^/^ ~ 3.7 times than for a = \ while the values 
of 5 = 1 are comparable to those for lower values of a. 

The liquid-vapor transition is of particular interest as the metastable liquid phase is thought 
to play a key role in the process of precipitation of solid protein crystals from solution.-i*2j^ In 
classical nucleation theory (CNT) a liquid droplet is treated in the capillary approximation so that 
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a droplet of radius R has excess free energy 

A^^—R^(o{pi)-0}{py))+47rR^r (12) 

where tt)(p) = /(p) — fip, f{p) is the bulk Helmholtz free energy per unit volume of the fluid at 
density p, /i is the chemical potential, 7 is the liquid- vapor surface tension at coexistence and p/ 
and py are the densities of the coexisting liquid and vapor respectively. In this model, the excess 
number of molecules in the droplet is 

An^—R\pi-p„). (13) 

The free energy has a maximum at the critical radius, Rc = ly/Ao), with Aco = 0){py) — co{pi), 
and a maximum, defining the barrier for nucleation, of AQ.max = Thus, from these simple 

considerations and the previously noted trends in the surface tension, we conclude that if all of 
the energy scales are proportional to Tc, then the nucleation barrier will decrease with increasing 
molecular size (5) and will vary weakly with the range of the potential, a. 

To investigate the liquid-vapor transition in more detail, we have determined the nucleation 
pathway (specifically, the minimum free energy path or MFEP) for nucleation of liquid droplets 
from the DFT to compare to the predictions of CNT for T = 0.4e or T/Tc = 0.84. The method of 
calculation is the same as described in detail in Refs.— >^ The goal is to identify a path between the 
initial state (pure vapor) and the final state (pure liquid) which is minimal with respect to variations 
perpendicular to the path (for a detailed explanation, see e.g.—). [figure] [6] []6 shows the free 
energy barrier and the size of the critical cluster as a function of supersaturation, S = ^'p"~^ where 

"coex 

the chemical potential at coexistence is j8 /Jcoex = —2.354 and the densities of the coexisting phases 
are Pv = 0.109 and p/ = 0.626. (Note that although the actual control variable in our calculations 
in the grand canonical ensemble is chemical potential, we use the more familiar definition of 



supersaturation in terms of the pressure.) [figure] [6] []6 shows that classical nucleation theory is in 



agreement with DFT for small super saturations but for large supersaturation, the nucleation barrier 



16 




Figure 6: (Color online) Panel (a) shows the free energy barrier for nucleation of dense liquid 
from the low-density phase as a function of supersaturation. Panel (b) shows the size of the critical 
cluster. In both cases, the full lines are from the DFT calculations and the broken lines are the 
prediction of CNT. 

calculated with DFT becomes very small, approaching zero as the vapor density approaches the 
spinodal, whereas CNT predicts a finite nucleation barrier at all densities. The figure also shows 
that the excess number of molecules in the critical cluster is also well described by CNT except 



at large supersaturations. [figure] [7] [] 7 showing the critical clusters at different supersaturations, 
illustrates the reason for the failure of CNT at large supersaturation. It shows that large clusters are 
indeed well described with the capillary model, having very narrow interfaces. However, at large 
supersaturation, the critical clusters are very small with most molecules affected by the interface 
and with central densities far below that of the bulk liquid. As a consequence, CNT fails to capture 



the very small barriers for nucleation. This final point is illustrated in [figure] [8] [] 8 which shows 
the nucleation pathway, i.e. the excess free energy as a function of cluster size, for three different 
supersaturations. At the smallest supersaturation, the barrier is well described by CNT. However, 
at large supersaturation, the shape of the nucleation barrier varies markedly from that assumed in 
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Figure 7: (Color online) The density distribution in the critical clusters for various values of the 
supersaturation. 

CNT. In particular, there is virtually no free energy penalty for the formation of small clusters. 

Given that we have access to the nucleation pathway, and not just the critical cluster, it is 
possible to directly evaluate the nucleation rate. Under the assumption of stationary nucleation and 
of treating the number of molecules as a continuous variable, the nucleation rate is given by the 
otherwise exact expression 

dn ^ ^ 



^-U/Wcwj '''' 

where f{n) is the monomer attachment rate for a cluster of size n and C{n) = CQexp{—Q.{n) /ksT) 
is the equilibrium distribution of cluster sizes.—""— The monomer attachment rate is assumed to 
be proportional to the surface area of a cluster and to the gas pressure, f{n) = fo{P/kBT)n^/^ with 
fo = 7'^'%^^ iji^^y^^ where 7 is the sticking probability for a monomer that collides with a cluster 
of size n (assumed here to be independent of n), c is a geometric factor and mq and vq are the mass 
and volume of a molecule so that the area of the cluster is cvq'^^h^/^.— [figure] [9] []9 shows the 



calculated nucleation rate as a function of the supersaturation. The nucleation rate estimated from 
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Figure 8: The excess free energy as a function of cluster size at three different values of supersatu- 
ration. The broken lines are a fit to the CNT expression Aa = -^AN+j^AN^/^ where the 

starred-quantities refer to the critical nucleus. 
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Figure 9: (Color online) Dimensionless nucleation rate as a function of the supersaturation. The 
scaled nucleation rate is J* = Ja^ / /q. 

DFT is significantly higher than that estimated from CNT, partly due to the exponential amplifica- 



tion of the differences in barrier height seen in [figure] [6] [] 6 and partly due to the differences in 
shape of the barrier. 



[figure] [6] [] 6 and [figure] [9] []9 both show non-monotonic behavior near the spinodal. In the 
later case, the implication is that there is an optimal supersaturation that gives the maximal nu- 
cleation rate. However, this is likely to be an artifact of the calculations. As discussed in detail 
by Wilemski and Li,— mean field models such as that used here, and simpler models such as the 
square gradient model, predict a divergent critical nucleus size at the spinodal due to the the unre- 
alistic mean-field equation of state. In any case, since the nucleation barrier is on the order of fc^T 
for S ~ 0.4, it is likely that the process of nucleation at higher supersaturations shares features of 
both nucleation and spinodal decomposition sometimes termed "spinodal nucleation"— in which 
case the nucleation characteristics (rate and critical nucleus) would not be described by the simple 
theory given here. There are, therefore, two issues affecting the interpretation of these results. The 
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first is tlie validity of mean field theory, which definitely breaks down near the spinodal-^^^^ and 
the second is that even if the maximum occurred in the region of validity of the theory, it would be 
masked by density fluctuations that would cause small volumes to become unstable with respect to 
spinodal decomposition. 



Conclusions 

In this paper, we have constructed a pair potential that allows us to move smoothly from a simple 
fluid (i.e. Lennard- Jones interaction) to the short-ranged, "hard Lennard- Jones" used as a model for 
globular protein interactions. The potential depends on two parameters: the size of the molecule, 
characterized by the excluded volume parameter 5, and the range of the potential, controlled by 
the dimensionless parameter a. We compared the liquid-vapor surface tension as determined from 
Monte Carlo simulations to DFT calculations for the case a = 1 , 5 and molecular radius ranging 
from zero to one (in Lennard- Jones units). Our results indicate that the combination of thermody- 
namic perturbation theory and the MC-VDW DFT model give a good quantitative description of 
the liquid- vapor equation of state and surface tension over a wide range of temperatures. Signifi- 
cant differences do appear as expected near the critical point since the mean-field equation of state 
is not accurate in this region. Our primary conclusion from these calculations is that it is the in- 
crease in excluded volume rather than the decrease in range of the potential that causes a dramatic 



decrease in the surface tension. In fact, from [figure] [4] []4 and [figure] [5] []5 , one sees that while 



the surface tension for zero hard-core radius is strongly affected by the range of the potential, the 
surface tension at the largest excluded volume is relatively insensitive to the range. 

For fixed 5, changing a changes the range of the potential. However, the model potential can 
in fact be rewritten in terms of a single dimensioness parameter, ad^ so that in absence of a cutoff, 
changing a and changing 5 are in some sense equivalent. The actual physical relevance of the 
two parameters is that 5 controls the size of the exluded volume while a controls the strength 
of the attractive interaction. Hence, one way to phrase our results is that changing the excluded 



21 



volume has much more effect on physical properties, such as surface tension, than does changing 
the strength of the attraction. These results suggest that increasing the excluded volume of a 
molecule leads to a decrease in surface tension and hence of the barrier to nucleation of the dense 
phase. Conformational changes can in fact be affected in some cases by changes in pH~""— and by 
light 3^ thus suggesting a means for taking advantage of this phenomena to control crystallization 
rates in such systems. 

For the case a = 50, we were not able to simulate the liquid- vapor interface due to the strong 
instability of the system to crystallization. Our theoretical calculations indicate similar behavior 
in the surface tension as found for smaller values of a. We have also calculated the nucleation 
pathway and thereby determined the barrier and rate of nucleation of liquid droplets from the 
vapor. We found quite low barriers, less than lOO^sT even at relatively low super saturation. Away 
from the spinodal. Classical Nucleation Theory gives a good description of the barrier height and 
the size of the critical nucleus. At higher supersaturations, as the spinodal is approached, the 
nucleation pathway calculated from DFT differs significantly from that predicted by CNT. DPT 
also predicts a vanishing nucleation barrier at the spinodal and non-monotonic behavior of the size 
of the critical nucleus and the nucleation rate near the spinodal. The shape of the barriers suggests 
that the excess free energy of small, sub-critical clusters is very small indicating that they might 
have relatively long lifetimes. However, the details of this picture become increasingly uncertain 
near the spinodal due to the limitations of mean field theory and in any case may not be observable 
as they occur in the region of "spinodal nucleation". 
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